#!/usr/bin/python

import sys, math, random, time
import numpy as np
from profilehooks import profile

# integrate with space saving
# change iteration to work with DCSC, DCSR
# implement MPI / mmap
# exploit symmetry: upper/lower triangular
#

def generate_matrix(n, r) :
    global A,B, C_mat, p_dist, C, h1, h2
    C, h1,h2 = {}, {},{}
    p_dist = np.zeros((p,), dtype=np.int)
#    A = B = np.load('stoc_mat_1000.npy')
    A = np.random.randint(r, size=(n, n))
    B = np.random.randint(r, size=(n, n))
    C_mat = [[0 for i in range(n)] for x in range(n)]

def val(h, x, hi):
    if not(x in h):
	random.seed(x + hi) # http://docs.python.org/library/random.htm. Hashvalue of x is determined by the random.seed
	h[x] = int(random.random() * p)
    return h[x]

def do_outer_product(h1val, h2val, h1_counts) :
    for pi in range(p) :
        for b_item in h2val :
            a_idx = (p + pi - b_item[0]) % p
            a_range = h1val[h1_counts[a_idx]:h1_counts[a_idx + 1]]

            p_dist[pi] += len(a_range)
            for a_item in a_range :
                key = (a_item[1][0], b_item[1][0])
                if key not in C :
                    C[key] = 0
                C[key] += b_item[1][1] * a_item[1][1]
def CRoP() :
    for a_k, col in enumerate(A.T) :
	h1_counts = [0 for x in range(p)]
	h1val, h2val = [], [] # list of tuples: (hashvalue of x, x)
	for a_i, item in enumerate(col) :
	    a_val = int(item)
	    hash1 = val(h1, a_i, 0)
	    h1_counts[hash1] += 1
	    h1val.append((hash1, (a_i, a_val)))

	for b_j, value in enumerate(B[a_k]) :
	    h2val.append((val(h2, b_j, 1), (b_j, int(value))))

	h1_counts = list(np.cumsum(h1_counts))
	h1_counts.insert(0, 0)

	h1val.sort()
	do_outer_product(h1val, h2val, h1_counts)

global p
p = 10
generate_matrix(100, 10)
CRoP()

for key, _val in C.iteritems():
    try :
	C_mat[key[0]][key[1]] = val
    except Exception as e :
        print e
	print key, _val

print "p_dist: ", p_dist
print "avg: ", np.average(p_dist)
print "mean: ", np.mean(p_dist)
print "std: ", np.std(p_dist)

p = 1
p_dist = np.zeros((p,), dtype=np.int)
C.clear()
h1.clear()
h2.clear()
CRoP()

print "p_dist: ", p_dist
print "avg: ", np.average(p_dist)
print "mean: ", np.mean(p_dist)
print "std: ", np.std(p_dist)

A = np.mat(A)
B = np.mat(B)
COUT = A * B
C_mat = np.mat(C_mat)

print (np.array_equal(COUT, C_mat))